home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / cagd_lib / cagdpris.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  14.7 KB  |  453 lines
  1. /******************************************************************************
  2. * CagdPrisa.c - piecewise ruled srf approximation and layout (prisa).          *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Apr. 93.                          *
  5. ******************************************************************************/
  6.  
  7. #include "cagd_loc.h"
  8.  
  9. static CagdSrfStruct *CagdPiecewiseRuledSrfAux(CagdSrfStruct *Srf,
  10.                            CagdBType ConsistentDir,
  11.                            CagdRType Epsilon,
  12.                            CagdSrfDirType Dir);
  13. static void CagdInterCircCirc(CagdRType Center1[3], CagdRType Radius1,
  14.                   CagdRType Center2[3], CagdRType Radius2,
  15.                   CagdRType Inter1[3], CagdRType Inter2[3]);
  16.  
  17. /******************************************************************************
  18. * Combine all the layout/prisa routines. Computes a piecewise ruled surface   *
  19. * approximation to a given set of surfaces with given Epsilon, and place them *
  20. * nicely on the plane, by approximating each ruled curve by SamplesPerCurve   *
  21. * samples. Dir controls the direction of ruled approximation, SpaceScale and  *
  22. * Offset controls the placement of the different planar pieces.              *
  23. ******************************************************************************/
  24. CagdSrfStruct *CagdAllPrisaSrfs(CagdSrfStruct *Srfs, int SamplesPerCurve,
  25.                     CagdRType Epsilon, CagdSrfDirType Dir,
  26.                     CagdVType Space)
  27. {
  28.     int SrfIndex = 1;
  29.     CagdSrfStruct *Srf,
  30.     *PrisaSrfsAll = NULL;
  31.     CagdVType Offset;
  32.  
  33.     for (Srf = Srfs; Srf != NULL; Srf = Srf -> Pnext, SrfIndex++) {
  34.     CagdSrfStruct *TSrf, *RSrf, *RuledSrfs;
  35.  
  36.     if (Epsilon > 0) {
  37.         RuledSrfs = CagdPiecewiseRuledSrfApprox(Srf, FALSE, Epsilon, Dir);
  38.  
  39.         Offset[0] = SrfIndex * Space[0];
  40.         Offset[1] = 0.0;
  41.         Offset[2] = 0.0;
  42.  
  43.         for (RSrf = RuledSrfs; RSrf != NULL; RSrf = RSrf -> Pnext) {
  44.         TSrf = CagdPrisaRuledSrf(RSrf, SamplesPerCurve,
  45.                      Space[1], Offset);
  46.         TSrf -> Pnext = PrisaSrfsAll;
  47.         PrisaSrfsAll = TSrf;
  48.         }
  49.  
  50.         CagdSrfFreeList(RuledSrfs);
  51.     }
  52.     else {
  53.         /* Return the 3D ruled surface approximation. */
  54.         RuledSrfs = CagdPiecewiseRuledSrfApprox(Srf, FALSE, -Epsilon, Dir);
  55.  
  56.         for (TSrf = RuledSrfs;
  57.          TSrf -> Pnext != NULL;
  58.          TSrf = TSrf -> Pnext);
  59.         TSrf -> Pnext = PrisaSrfsAll;
  60.         PrisaSrfsAll = RuledSrfs;
  61.     }
  62.     }
  63.  
  64.     return PrisaSrfsAll;
  65. }
  66.  
  67. /******************************************************************************
  68. * Constructs a piecewise ruled surface approximation to the given surface in  *
  69. * the given direction, that is close to the surface to within Epsilon.          *
  70. *   If ConsitentDir then ruled surface parametrization is set to be the same  *
  71. * same as original surface. Otherwise, ruling dir is always CAGD_CONST_V_DIR. *
  72. ******************************************************************************/
  73. CagdSrfStruct *CagdPiecewiseRuledSrfApprox(CagdSrfStruct *Srf,
  74.                        CagdBType ConsistentDir,
  75.                        CagdRType Epsilon,
  76.                        CagdSrfDirType Dir)
  77. {
  78.     CagdSrfStruct *RuledSrfs;
  79.  
  80.     switch (Srf -> PType) {
  81.     case CAGD_PT_E3_TYPE:
  82.     case CAGD_PT_P3_TYPE:
  83.         Srf = CagdSrfCopy(Srf);
  84.         break;
  85.     case CAGD_PT_P1_TYPE:
  86.     case CAGD_PT_P2_TYPE:
  87.     case CAGD_PT_P4_TYPE:
  88.     case CAGD_PT_P5_TYPE:
  89.         Srf = CagdCoerceSrfTo(Srf, CAGD_PT_P3_TYPE);
  90.         break;
  91.     case CAGD_PT_E1_TYPE:
  92.     case CAGD_PT_E2_TYPE:
  93.     case CAGD_PT_E4_TYPE:
  94.     case CAGD_PT_E5_TYPE:
  95.         Srf = CagdCoerceSrfTo(Srf, CAGD_PT_E3_TYPE);
  96.         break;
  97.     default:
  98.         FATAL_ERROR(CAGD_ERR_UNSUPPORT_PT);
  99.         Srf = NULL;
  100.         break;
  101.     }
  102.  
  103.     RuledSrfs = CagdPiecewiseRuledSrfAux(Srf, ConsistentDir, Epsilon, Dir);
  104.     CagdSrfFree(Srf);
  105.     return RuledSrfs;
  106. }
  107.  
  108. /******************************************************************************
  109. * Constructs a piecewise ruled surface approximation to the given surface in  *
  110. * the given direction, that is close to the surface to within Epsilon.          *
  111. * Surface is assumed to have point types E3 or P3 only.                  *
  112. ******************************************************************************/
  113. static CagdSrfStruct *CagdPiecewiseRuledSrfAux(CagdSrfStruct *Srf,
  114.                            CagdBType ConsistentDir,
  115.                            CagdRType Epsilon,
  116.                            CagdSrfDirType Dir)
  117. {
  118.     CagdSrfStruct *DiffSrf, *DistSqrSrf,
  119.     *RuledSrf = CagdSrfCopy(Srf);
  120.     CagdRType *XPts, *WPts, UMin, UMax, VMin, VMax,
  121.     **Points = RuledSrf -> Points,
  122.     MaxError = 0.0,
  123.     t = 0.0;
  124.     int i, j, k, Index,
  125.     ULength = RuledSrf -> ULength,
  126.     VLength = RuledSrf -> VLength;
  127.     PointType E3PtStart, E3PtEnd, E3Pt;
  128.  
  129.     switch (Dir) {
  130.     case CAGD_CONST_U_DIR:
  131.         for (j = 0; j < VLength; j++) {
  132.         /* First order approximation to the ratios of the   */
  133.         /* distance of interior point to the end points.    */
  134.         CagdCoerceToE3(E3PtStart, Points,
  135.                    CAGD_MESH_UV(RuledSrf, ULength / 2, 0),
  136.                    Srf -> PType);
  137.         CagdCoerceToE3(E3PtEnd, Points,
  138.                    CAGD_MESH_UV(RuledSrf, ULength / 2,
  139.                                   VLength - 1),
  140.                    Srf -> PType);
  141.         CagdCoerceToE3(E3Pt, Points,
  142.                    CAGD_MESH_UV(RuledSrf, ULength / 2, j),
  143.                    Srf -> PType);
  144.         PT_SUB(E3PtStart, E3PtStart, E3PtEnd);
  145.         PT_SUB(E3Pt, E3Pt, E3PtEnd);
  146.         t = PT_LENGTH(E3PtStart);
  147.         t = APX_EQ(t, 0.0) ? 0.5 : PT_LENGTH(E3Pt) / t;
  148.  
  149.         for (i = 0; i < ULength; i++) {
  150.             CagdCoerceToE3(E3PtStart, Points,
  151.                    CAGD_MESH_UV(RuledSrf, i, 0),
  152.                    Srf -> PType);
  153.             CagdCoerceToE3(E3PtEnd, Points,
  154.                    CAGD_MESH_UV(RuledSrf, i, VLength - 1),
  155.                    Srf -> PType);
  156.             PT_BLEND(E3Pt, E3PtStart, E3PtEnd, t);
  157.  
  158.             Index = CAGD_MESH_UV(RuledSrf, i, j);
  159.             if (CAGD_IS_RATIONAL_PT(RuledSrf -> PType)) {
  160.             for (k = 0; k < 3; k++)
  161.                 Points[k + 1][Index] = 
  162.                 E3Pt[k] * Points[0][Index];
  163.             }
  164.             else {
  165.             for (k = 0; k < 3; k++)
  166.                 Points[k + 1][Index] = E3Pt[k];
  167.             }
  168.         }
  169.         }
  170.         break;
  171.     case CAGD_CONST_V_DIR:
  172.         for (i = 0; i < ULength; i++) {
  173.         /* First order approximation to the ratios of the   */
  174.         /* distance of interior point to the end points.    */
  175.         CagdCoerceToE3(E3PtStart, Points,
  176.                    CAGD_MESH_UV(RuledSrf, 0, VLength / 2),
  177.                    Srf -> PType);
  178.         CagdCoerceToE3(E3PtEnd, Points,
  179.                    CAGD_MESH_UV(RuledSrf, ULength - 1,
  180.                               VLength / 2),
  181.                    Srf -> PType);
  182.         CagdCoerceToE3(E3Pt, Points,
  183.                    CAGD_MESH_UV(RuledSrf, i, VLength / 2),
  184.                    Srf -> PType);
  185.         PT_SUB(E3PtStart, E3PtStart, E3PtEnd);
  186.         PT_SUB(E3Pt, E3Pt, E3PtEnd);
  187.         t = PT_LENGTH(E3PtStart);
  188.         t = APX_EQ(t, 0.0) ? 0.5 : PT_LENGTH(E3Pt) / t;
  189.  
  190.         for (j = 0; j < VLength; j++) {
  191.             CagdCoerceToE3(E3PtStart, Points,
  192.                    CAGD_MESH_UV(RuledSrf, 0, j),
  193.                    Srf -> PType);
  194.             CagdCoerceToE3(E3PtEnd, Points,
  195.                    CAGD_MESH_UV(RuledSrf, ULength - 1, j),
  196.                    Srf -> PType);
  197.             PT_BLEND(E3Pt, E3PtStart, E3PtEnd, t);
  198.  
  199.             Index = CAGD_MESH_UV(RuledSrf, i, j);
  200.             if (CAGD_IS_RATIONAL_PT(RuledSrf -> PType)) {
  201.             for (k = 0; k < 3; k++)
  202.                 Points[k + 1][Index] = 
  203.                 E3Pt[k] * Points[0][Index];
  204.             }
  205.             else {
  206.             for (k = 0; k < 3; k++)
  207.                 Points[k + 1][Index] = E3Pt[k];
  208.             }
  209.         }
  210.         }
  211.         break;
  212.     default:
  213.         FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  214.         break;
  215.     }
  216.  
  217.     DiffSrf = CagdSrfSub(Srf, RuledSrf);
  218.     CagdSrfFree(RuledSrf);
  219.     DistSqrSrf = CagdSrfDotProd(DiffSrf, DiffSrf);
  220.     CagdSrfFree(DiffSrf);
  221.     XPts = DistSqrSrf -> Points[1];
  222.     WPts = CAGD_IS_RATIONAL_PT(DistSqrSrf -> PType) ?
  223.                     DistSqrSrf -> Points[0] : NULL;
  224.  
  225.     for (i = DistSqrSrf -> ULength * DistSqrSrf -> VLength; i > 0; i--) {
  226.     CagdRType
  227.         V = WPts != NULL ? *XPts++ / *WPts++ : *XPts++;
  228.  
  229.     if (MaxError < V)
  230.         MaxError = V;
  231.     }
  232.     CagdSrfFree(DistSqrSrf);
  233.  
  234.     if (MaxError > SQR(Epsilon)) {
  235.     /* Subdivide and try again. */
  236.     CagdSrfStruct *Srf1, *Srf2, *RuledSrf1, *RuledSrf2;
  237.  
  238.     CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  239.     t = Dir == CAGD_CONST_V_DIR ? (UMax + UMin) / 2 : (VMax + VMin) / 2;
  240.     Srf1 = CagdSrfSubdivAtParam(Srf, t, CAGD_OTHER_DIR(Dir));
  241.     Srf2 = Srf1 -> Pnext;
  242.     Srf1 -> Pnext = NULL;
  243.  
  244.     RuledSrf1 = CagdPiecewiseRuledSrfAux(Srf1, ConsistentDir, Epsilon, Dir);
  245.     RuledSrf2 = CagdPiecewiseRuledSrfAux(Srf2, ConsistentDir, Epsilon, Dir);
  246.     CagdSrfFree(Srf1);
  247.     CagdSrfFree(Srf2);
  248.  
  249.     for (Srf1 = RuledSrf1; Srf1 -> Pnext != NULL; Srf1 = Srf1 -> Pnext);
  250.     Srf1 -> Pnext = RuledSrf2;
  251.     return RuledSrf1;
  252.     }
  253.     else {
  254.     /* Returns the ruled surface as a linear in the ruled direction. */
  255.     CagdCrvStruct
  256.         *Crv1 = CagdCrvFromMesh(Srf, 0, CAGD_OTHER_DIR(Dir)),
  257.         *Crv2 = CagdCrvFromMesh(Srf,
  258.                     Dir == CAGD_CONST_V_DIR ? ULength - 1
  259.                                 : VLength - 1,
  260.                     CAGD_OTHER_DIR(Dir));
  261.  
  262.     RuledSrf = CagdRuledSrf(Crv1, Crv2, 2, 2);
  263.     if (ConsistentDir && Dir == CAGD_CONST_V_DIR) {
  264.         /* Needs to reverse the ruled surface so it matches Srf. */
  265.         CagdSrfStruct
  266.         *TSrf = CagdSrfReverse2(RuledSrf);
  267.  
  268.         CagdSrfFree(RuledSrf);
  269.         RuledSrf = TSrf;
  270.     }
  271.  
  272.     if (CAGD_IS_BSPLINE_SRF(Srf)) {
  273.         /* Updates the knot vector domain. */
  274.         CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  275.         if (Dir == CAGD_CONST_V_DIR)
  276.             BspKnotAffineTrans(RuledSrf -> UKnotVector,
  277.                    RuledSrf -> ULength + RuledSrf -> UOrder,
  278.                    UMin, UMax - UMin);
  279.         else
  280.             BspKnotAffineTrans(RuledSrf -> VKnotVector,
  281.                    RuledSrf -> VLength + RuledSrf -> VOrder,
  282.                    VMin, VMax - VMin);
  283.     }
  284.     CagdCrvFree(Crv1);
  285.     CagdCrvFree(Crv2);
  286.     return RuledSrf;
  287.     }
  288. }
  289.  
  290. /******************************************************************************
  291. * Layout a ruled surface, by approximating it as a set of polygons.          *
  292. *   The given ruled surface might be non-developable, in which case          *
  293. * approximation will be of a surface with no twist.                  *
  294. *   The ruled surface is assumed to be constructed using CagdRuledSrf and the *
  295. * ruled direction is always CAGD_CONST_V_DIR.                      *
  296. ******************************************************************************/
  297. CagdSrfStruct *CagdPrisaRuledSrf(CagdSrfStruct *Srf, int SamplesPerCurve,
  298.                     CagdRType Space, CagdVType Offset)
  299. {
  300.     int i, j,
  301.     VLength = Srf -> VLength;
  302.     CagdRType Angle, PtTmp1[3], PtTmp2[3], PtTmp3[3];
  303.     CagdCrvStruct
  304.     *Crv1 = CagdCrvFromMesh(Srf, 0, CAGD_CONST_V_DIR),
  305.     *Crv2 = CagdCrvFromMesh(Srf, VLength - 1, CAGD_CONST_V_DIR);
  306.     CagdPolylineStruct
  307.     *Poly1 = CagdCrv2Polyline(Crv1, SamplesPerCurve),
  308.     *Poly2 = CagdCrv2Polyline(Crv2, SamplesPerCurve),
  309.     *Poly1Prisa = CagdPolylineNew(Poly1 -> Length),
  310.     *Poly2Prisa = CagdPolylineNew(Poly2 -> Length);
  311.     CagdPtStruct
  312.         *Pt1 = Poly1 -> Polyline,
  313.         *Pt2 = Poly2 -> Polyline,
  314.         *Pt1Prisa = Poly1Prisa -> Polyline,
  315.         *Pt2Prisa = Poly2Prisa -> Polyline;
  316.     CagdMType Mat1, Mat;
  317.     CagdBBoxStruct BBox;
  318.  
  319.     CagdCrvFree(Crv1);
  320.     CagdCrvFree(Crv2);
  321.  
  322.     /* Anchor the location of the first line. */
  323.     for (j = 0; j < 3; j++) {
  324.     Pt1Prisa -> Pt[j] = 0.0;
  325.     Pt2Prisa -> Pt[j] = 0.0;
  326.     }
  327.     PT_SUB(PtTmp1, Pt1 -> Pt, Pt2 -> Pt);
  328.     Pt2Prisa -> Pt[1] = PT_LENGTH(PtTmp1);
  329.  
  330.     /* Move alternately along Poly1 and Poly2 and anchor the next point of */
  331.     /* the next triangle using the distances to current Pt1 and Pt2.       */
  332.     for (i = 2; i < Poly1 -> Length + Poly2 -> Length; i++) {
  333.     CagdRType Dist1, Dist2, Inter1[3], Inter2[3],
  334.         *NextPt = i & 0x01 ? Pt1[1].Pt : Pt2[1].Pt;
  335.  
  336.     /* Compute distance from previous two locations to new location. */ 
  337.     PT_SUB(PtTmp1, Pt1 -> Pt, NextPt);
  338.     Dist1 = PT_LENGTH(PtTmp1);
  339.     PT_SUB(PtTmp1, Pt2 -> Pt, NextPt);
  340.     Dist2 = PT_LENGTH(PtTmp1);
  341.  
  342.     /* Find the (two) intersection points of circles with radii Dist?. */
  343.     CagdInterCircCirc(Pt1Prisa -> Pt, Dist1, Pt2Prisa -> Pt, Dist2,
  344.               Inter1, Inter2);
  345.  
  346.     /* Find which of the two intersection points is the "good" one. */
  347.     PT_SUB(PtTmp1, Inter1, Pt1Prisa -> Pt);
  348.     PT_SUB(PtTmp2, Inter1, Pt2Prisa -> Pt);
  349.     CROSS_PROD(PtTmp3, PtTmp1, PtTmp2);
  350.     if (i & 0x01) {
  351.         Pt1Prisa++;
  352.         for (j = 0; j < 3; j++)
  353.         Pt1Prisa -> Pt[j] = (PtTmp3[2] > 0 ? Inter2[j] : Inter1[j]);
  354.         Pt1++;
  355.     }
  356.     else {
  357.         Pt2Prisa++;
  358.         for (j = 0; j < 3; j++)
  359.         Pt2Prisa -> Pt[j] = (PtTmp3[2] > 0 ? Inter2[j] : Inter1[j]);
  360.         Pt2++;
  361.     }
  362.     }
  363.  
  364.     /* Save centering location so we can orient the resulting data nicely. */
  365.     PT_COPY(PtTmp1, Poly1Prisa -> Polyline[Poly1Prisa -> Length / 2].Pt);
  366.     PT_COPY(PtTmp2, Poly2Prisa -> Polyline[Poly2Prisa -> Length / 2].Pt);
  367.     PT_SUB(PtTmp3, PtTmp2, PtTmp1);
  368.  
  369.     Crv1 = CnvrtPolyline2LinBsplineCrv(Poly1Prisa);
  370.     Crv2 = CnvrtPolyline2LinBsplineCrv(Poly2Prisa);
  371.     CagdPolylineFree(Poly1);
  372.     CagdPolylineFree(Poly2);
  373.     CagdPolylineFree(Poly1Prisa);
  374.     CagdPolylineFree(Poly2Prisa);
  375.  
  376.     Srf = CagdRuledSrf(Crv1, Crv2, 2, 2);
  377.     CagdCrvFree(Crv1);
  378.     CagdCrvFree(Crv2);
  379.  
  380.     /* Translate PtTmp1 to the origin. */
  381.     MatGenMatTrans(-PtTmp1[0], -PtTmp1[1], -PtTmp1[2], Mat);
  382.  
  383.     /* Rotate PtTmp3 direction to the +Y direction. */
  384.     Angle = atan2(PtTmp3[1], PtTmp3[0]);
  385.     MatGenMatRotZ1(M_PI / 2 - Angle, Mat1);
  386.  
  387.     MatMultTwo4by4(Mat, Mat, Mat1);
  388.  
  389.     CagdSrfMatTransform(Srf, Mat);
  390.  
  391.     /* Translate by the Offset. */
  392.     CagdSrfBBox(Srf, &BBox);
  393.     MatGenMatTrans(Offset[0], Offset[1] - BBox.Min[1], Offset[2], Mat);
  394.     Offset[1] += (BBox.Max[1] - BBox.Min[1]) + Space;
  395.  
  396.     CagdSrfMatTransform(Srf, Mat);
  397.  
  398.     return Srf;
  399. }
  400.  
  401. /******************************************************************************
  402. * Find the two intersection points of the given two planar circles. It is     *
  403. * assumed intersection exists.                              *
  404. ******************************************************************************/
  405. static void CagdInterCircCirc(CagdRType Center1[3], CagdRType Radius1,
  406.                   CagdRType Center2[3], CagdRType Radius2,
  407.                   CagdRType Inter1[3], CagdRType Inter2[3])
  408. {
  409.     CagdRType A, B, C, Delta,
  410.     a = Center2[0] - Center1[0],
  411.     b = Center2[1] - Center1[1],
  412.     c = (SQR(Radius1) - SQR(Radius2) +
  413.          SQR(Center2[0]) - SQR(Center1[0]) +
  414.          SQR(Center2[1]) - SQR(Center1[1])) / 2;
  415.  
  416.     if (PT_APX_EQ(Center1, Center2)) {
  417.     Inter1[0] = Inter2[0] = Radius1;
  418.     Inter1[1] = Inter2[1] = 0.0;
  419.     }
  420.     else if (ABS(a) > ABS(b)) {
  421.     /* Solve for Y first. */
  422.     A = SQR(b) / SQR(a) + 1;
  423.     B = 2 * (b * Center1[0] / a - b * c / SQR(a) - Center1[1]);
  424.     C = SQR(c / a) + SQR(Center1[0]) + SQR(Center1[1])
  425.                 -2 * c * Center1[0] / a - SQR(Radius1);
  426.     Delta = SQR(B) - 4 * A * C;
  427.     if (Delta < 0)  /* If no solution, do something almost reasonable. */
  428.         Delta = 0;
  429.     Inter1[1] = (-B + sqrt(Delta)) / (2 * A);
  430.     Inter2[1] = (-B - sqrt(Delta)) / (2 * A);
  431.  
  432.     Inter1[0] = (c - b * Inter1[1]) / a;
  433.     Inter2[0] = (c - b * Inter2[1]) / a;
  434.     }
  435.     else {
  436.     /* Solve for X first. */
  437.     A = SQR(a) / SQR(b) + 1;
  438.     B = 2 * (a * Center1[1] / b - a * c / SQR(b) - Center1[0]);
  439.     C = SQR(c / b) + SQR(Center1[0]) + SQR(Center1[1])
  440.                 -2 * c * Center1[1] / b - SQR(Radius1);
  441.     Delta = SQR(B) - 4 * A * C;
  442.     if (Delta < 0)  /* If no solution, do something almost reasonable. */
  443.         Delta = 0;
  444.     Inter1[0] = (-B + sqrt(Delta)) / (2 * A);
  445.     Inter2[0] = (-B - sqrt(Delta)) / (2 * A);
  446.  
  447.     Inter1[1] = (c - a * Inter1[0]) / b;
  448.     Inter2[1] = (c - a * Inter2[0]) / b;
  449.     }
  450.  
  451.     Inter1[2] = Inter2[2] = 0.0;
  452. }
  453.